package scu.maqiang.fes;

import scu.maqiang.mesh.*;
import scu.maqiang.numeric.CG;
import scu.maqiang.numeric.DoubleIntegerPair;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.ParamCheck;
import scu.maqiang.numeric.SRMatrix;

public class FES3T41 extends FES {

	public FES3T41(Mesh3T4 mesh3D) {
		super(mesh3D, 1);
		eleStiff = FES3T41::elementStiff;
		eleMass = FES3T41::elementMass;
		eleSource = FES3T41::elementSource;
		eleFlux = FES3T41::elementHeatFlux;
	}

	public int getNdof() {
		return ndof;
	}

	public double getValue(double[] u, double[] xyz) {
		return getValue0(u, xyz, CG::pointInT4);
	}
	
	public double getValue(double[] u, double[] xyz, int[] eleIdx, double[] eleVal) {
		return getValue0(u, xyz, eleIdx, eleVal, CG::pointInT4);
//		ParamCheck.checkEqual(u.length, this.ndof);
//		ParamCheck.checkEqual(eleIdx.length, this.nPerEle);
//		ParamCheck.checkEqual(eleVal.length, this.nPerEle);
//		double[] uu = new double[4];
//		double value = 0;
//		double[][] coord = new double[4][3];
//		for (int i = 0, n = mesh.getNt(); i < n; i++) {
//			int[] ele = mesh.getElement(i);
//			mesh.getCoordPerElement(ele, coord);
//			CG.T4ShapeFunction(coord, xyz, eleVal);
//			if (Math.abs(MVO.sum(eleVal) - 1.0) < 1.0e-10) {
//				for (int j = 0; j < 4; j++) {
//					uu[j] = u[ele[j]];
//					eleIdx[j] = ele[j];
//					value += uu[j] * eleVal[j];
//				}
//				return value;
//			}
//		}
//		throw new IllegalArgumentException("The Point ( " + Arrays.toString(xyz) + " ) is not in the domain!");
	}

	public void getValue(double[] u, double[] xyz, DoubleIntegerPair data, int[] eleIdx, double[] eleVal) {
		getValue0(u, xyz, data, eleIdx, eleVal, CG::pointInT4);
//		ParamCheck.checkEqual(u.length, this.ndof);
//		ParamCheck.checkEqual(eleIdx.length, this.nPerEle);
//		ParamCheck.checkEqual(eleVal.length, this.nPerEle);
//		double value = 0;
//		double[][] coord = new double[4][3];
//		int[] ele = mesh.getElement(data.second);
//		mesh.getCoordPerElement(ele, coord);
//		CG.T4ShapeFunction(coord, xyz, eleVal);
//		if (Math.abs(MVO.sum(eleVal) - 1.0) < 1.0e-10) {
//			for (int j = 0; j < 4; j++) {
//				// eleIdx[j] = ele.idx[j];
//				value += u[ele[j]] * eleVal[j];
//			}
//			data.first = value;
//			return;
//		} else {
//			for (int i = 0, n = mesh.getNt(); i < n; i++) {
//				ele = mesh.getElement(i);
//				mesh.getCoordPerElement(ele, coord);
//				CG.T4ShapeFunction(coord, xyz, eleVal);
//				if (Math.abs(MVO.sum(eleVal) - 1.0) < 1.0e-10) {
//					for (int j = 0; j < 4; j++) {
//						eleIdx[j] = ele[j];
//						value += u[ele[j]] * eleVal[j];
//					}
//					data.first = value;
//					data.second = i;
//					return;
//				}
//			}
//
//		}
//		throw new IllegalArgumentException("The Point ( " + Arrays.toString(xyz) + " ) is not in the domain!");
	}

	public static void elementStiff(int[] ele, double[][] xyz, double[][] kxx, BVPType tp, double[][] eleMatrix) {
		double[][] grad = new double[4][3];
//		System.out.println(MVO.toString(xyz));
		double volume = CG.gradient3DT4(xyz, grad);
//		System.out.println(volume);
//		System.out.println(MVO.toString(grad));
		double Af;
		double[][] temp;
		switch (tp) {
		case COMMON:
			Af = volume * MVO.average(kxx[0]);
			for (int i = 0; i < 4; i++) {
				for (int j = 0; j < 4; j++) {
					eleMatrix[i][j] = (grad[i][0] * grad[j][0] + grad[i][1] * grad[j][1] + grad[i][2] * grad[j][2])
							* Af;
				}
			}
			break;
		case CONSTITUTE_MATRIX_COMMON:
			temp = MVO.matmul(grad, kxx);
			for (int i = 0; i < 4; i++) {
				for (int j = 0; j < 4; j++) {
					eleMatrix[i][j] = 0.0;
					for (int k = 0; k < 3; k++) {
						eleMatrix[i][j] += temp[i][k] * grad[j][k] * volume;
					}
				}
			}
			break;
		default:
			throw new IllegalArgumentException("Wrong BVPType!");
		}
//		System.out.println(MVO.toString(eleMatrix));
//		System.exit(0);
	}

	public static void elementMass(int[] ele, double[][] xyz, double[][] coef, BVPType tp, double[][] elemMatrix) {
		double volume = CG.volume3DT4(xyz);
		double[] ccoef = coef[0];
		double f1 = ccoef[0];
		double f2 = ccoef[1];
		double f3 = ccoef[2];
		double f4 = ccoef[3];
		System.out.println(f1 + "\t" + f2 + "\t"+ volume);
		if (tp == BVPType.CONSIST_MASS) {
			elemMatrix[0][0] = (3 * f1 + f2 + f3 + f4) * volume / 60;
			elemMatrix[0][1] = elemMatrix[1][0] = (2 * f1 + 2 * f2 + f3 + f4) * volume / 120;
			elemMatrix[0][2] = elemMatrix[2][0] = (2 * f1 + f2 + 2 * f3 + f4) * volume / 120;
			elemMatrix[0][3] = elemMatrix[3][0] = (2 * f1 + f2 + f3 + 2 * f4) * volume / 120;
			elemMatrix[1][1] = (f1 + 3 * f2 + f3 + f4) * volume / 60;
			elemMatrix[1][2] = elemMatrix[2][1] = (f1 + 2 * f2 + 2 * f3 + f4) * volume / 120;
			elemMatrix[1][3] = elemMatrix[3][1] = (f1 + 2 * f2 + f3 + 2 * f4) * volume / 120;
			elemMatrix[2][2] = (f1 + f2 + 3 * f3 + f4) * volume / 60;
			elemMatrix[2][3] = elemMatrix[3][2] = (f1 + f2 + 2 * f3 + 2 * f4) * volume / 120;
			elemMatrix[3][3] = (f1 + f2 + f3 + 3 * f4) * volume / 60;
		} else if (tp == BVPType.LUMP_MASS) {
			elemMatrix[0][0] = (2 * f1 + f2 + f3 + f4) * volume / 20;
			elemMatrix[1][1] = (f1 + 2 * f2 + f3 + f4) * volume / 20;
			elemMatrix[2][2] = (f1 + f2 + 2 * f3 + f4) * volume / 20;
			elemMatrix[3][3] = (f1 + f2 + f3 + 2 * f4) * volume / 20;
		}
//		System.out.println(MVO.toString(elemMatrix));
//		System.exit(0);
	}
	
	public void assembleHeatThetaScheme(ScalarFunc[] coefFunc, double[] param, BVPType[] tp, SRMatrix LMatrix, SRMatrix RMatrix) {
		double[][] kxx = new double[1][nPerEle];
		double[][] rho = new double[1][nPerEle];
		double[][] eleMatrixStiff = new double[dofPerElement][dofPerElement];
		double[][] eleMatrixMass = new double[dofPerElement][dofPerElement];
		double[][] eleLMatrix = new double[dofPerElement][dofPerElement];
		double[][] eleRMatrix = new double[dofPerElement][dofPerElement];
		double dt = param[0];
		double theta = param[1];
		for (int i = 0; i < nE; i++) {
			int[] ele = mesh.getElement(i);
			mesh.getCoordPerElement(ele, coordRegion);
			int label = mesh.getElementLabel(i);
			for (int j = 0, nj = nPerEle; j < nj; j++) {
				kxx[0][j] = coefFunc[0].action(coordRegion[j], label, param);
				rho[0][j] = coefFunc[1].action(coordRegion[j], label, param);
			}
			elementStiff(ele, coordRegion, kxx, tp[0], eleMatrixStiff);
			elementMass(ele, coordRegion, rho, tp[1], eleMatrixMass);
			for(int ii = 0; ii < dofPerElement; ii++ ) {
				for(int jj = 0; jj < dofPerElement; jj++) {
					eleLMatrix[ii][jj] = eleMatrixMass[ii][jj] + theta * dt * eleMatrixStiff[ii][jj];
					eleRMatrix[ii][jj] = eleMatrixMass[ii][jj] - (1 - theta) * dt * eleMatrixStiff[ii][jj];
				}
			}
			LMatrix.assemble(ele, eleLMatrix);
			RMatrix.assemble(ele, eleRMatrix);
		}
	}
	public void assembleHeatThetaScheme(double[] coefFunc, BVPType[] tp, SRMatrix LMatrix, SRMatrix RMatrix) {
		assembleHeatThetaScheme0(coefFunc, tp, FES3T41::elementStiff, FES3T41::elementMass, LMatrix, RMatrix);
	}


	public static void elementSource(int[] ele, double[][] xyz, double[][] f, BVPType tp, double[] eleF) {
		double area = CG.volume3DT4(xyz);
		double sumF = MVO.sum(f[0]);
		switch (tp) {
		case COMMON:
			for (int i = 0; i < 4; i++) {
				eleF[i] = (f[0][i] + sumF) * area / 20.0;
			}
			break;
		default:
			throw new IllegalArgumentException("Wrong BVPType!");
		}

	}
	
	
	public void assembleHeatConvectSurface(ScalarFunc[] coefFunc, double[] param, BVPType tp, SRMatrix A, int... bdLabel) {
		assembleBoundaryMatrix(coefFunc, param, FES3T41::elementConvectSurface, tp, A, bdLabel);
	}
	
	public void assembleHeatConvectSurface(double[] constConvect, BVPType tp, SRMatrix A, int... bdLabel) {
		assembleBoundaryMatrix(constConvect, FES3T41::elementConvectSurface, tp, A, bdLabel);
	}
	
	public void assembleHeatConvectSurface(double[][] vectorConvect, BVPType tp, SRMatrix A, int... bdLabel) {
		assembleBoundaryMatrix(vectorConvect, FES3T41::elementConvectSurface, tp, A, bdLabel);
	}
	
	public static void elementConvectSurface(int[] ele, double[][] xyz, double[][] h, BVPType tp, double[][] eleMatrix) {
	    double area = CG.area3DT3(xyz);
	    double[] g = h[0];
	    if (tp == BVPType.CONSIST_MASS) {
	    	area /= 60.0;
	    	eleMatrix[0][0] = 2 * (3 * g[0] + g[1] + g[2]) * area;
	    	eleMatrix[1][1] = 2 * (g[0] + 3 * g[1] + g[2]) * area;
	    	eleMatrix[2][2] = 2 * (g[0] + g[1] + 3 * g[2]) * area;
	    	eleMatrix[1][0] = eleMatrix[0][1] = (2 * g[0] + 2 * g[1] + g[2]) * area;
	    	eleMatrix[2][0] = eleMatrix[0][2] = (2 * g[0] + g[1] + 2 * g[2]) * area;
	    	eleMatrix[2][1] = eleMatrix[1][2] = (g[0] + 2 * g[1] + 2 * g[2]) * area;
	    } else if (tp == BVPType.LUMP_MASS) {
	    	area /= 12.0;
	    	eleMatrix[0][0] = (2 * g[0] + g[1] + g[1]) * area;
	    	eleMatrix[1][1] = (g[0] + 2 * g[1] + g[2]) * area;
	    	eleMatrix[2][2] = (g[0] + g[1] + 2 * g[2]) * area;
	    } else
	    	throw new IllegalArgumentException("Wrong BVPType!");
	}

	public static void elementHeatFlux(int[] ele, double[][] coord, double[][] g, BVPType tp, double[] k) {
		
		double[] gg= g[0];
		double area = CG.area3DT3(coord);
		area /= 12.0;
		double sumG = MVO.sum(gg);
		if (tp == BVPType.COMMON) {
			k[0] = (sumG + gg[0]) * area;
			k[1] = (sumG + gg[1]) * area;
			k[2] = (sumG + gg[2]) * area;
		} else
			throw new IllegalArgumentException("Wrong BVPType!");
//	    System.out.println(area * 12.0);
//	    System.out.println(Arrays.toString(k));
	}
	
	public void assembleSourceDerivative(ScalarFunc[] coefFunc, double[] param, BVPType tp, double[][] RHS) {
		assembleGlobalVector(coefFunc, param, this::elementSourceDerivative, tp, RHS);
	}
	
	public void elementSourceDerivative(int[] ele, double[][] xyz, double[][] coef, BVPType tp, double[][] eleXY) {
		double[][] grad = new double[4][3];
		double volume = CG.gradient3DT4(xyz, grad);
		double temp = MVO.sum(coef[0]) * volume / 4.0;

		eleXY[0][0] = -temp * grad[0][0];
		eleXY[0][1] = -temp * grad[1][0];
		eleXY[0][2] = -temp * grad[2][0];
		eleXY[0][3] = -temp * grad[3][0];
		eleXY[1][0] = -temp * grad[0][1];
		eleXY[1][1] = -temp * grad[1][1];
		eleXY[1][2] = -temp * grad[2][1];
		eleXY[1][3] = -temp * grad[3][1];
		eleXY[2][0] = -temp * grad[0][2];
		eleXY[2][1] = -temp * grad[1][2];
		eleXY[2][2] = -temp * grad[2][2];
		eleXY[2][3] = -temp * grad[3][2];
	}
	
	public void computeGradient(double[] u, double[][] gradU) {
		int len = u.length;
		ParamCheck.checkEqual(len, this.ndof);
		double[][] coord = new double[4][3];
		double[][] grad = new double[4][3];
		double[] volCount = new double[len];
		double[] eleU = new double[4];
		MVO.fill(gradU, 0.0);
		for(int i = 0, n = mesh.getNt(); i < n; i++) {
			int[] ele = mesh.getElement(i);
			mesh.getCoordPerElement(ele, coord);
			for(int j = 0; j < 4; j++) {
				eleU[j] = u[ele[j]];
			}
			double volume = CG.gradient3DT4(coord, grad);
			double[] dx = MVO.col(grad, 0);
			double[] dy = MVO.col(grad, 1);
			double[] dz = MVO.col(grad, 2);
			double dxValue = MVO.dot_product(eleU, dx);
			double dyValue = MVO.dot_product(eleU, dy);
			double dzValue = MVO.dot_product(eleU, dz);
			for(int j = 0; j < 4; j++) {
				gradU[0][ele[j]] += dxValue * volume;
				gradU[1][ele[j]] += dyValue * volume;
				gradU[2][ele[j]] += dzValue * volume;
				volCount[ele[j]] += volume;
			}
		}
		
		for(int i = 0, n = mesh.getNv(); i < n; i++) {
			gradU[0][i] /= volCount[i];
			gradU[1][i] /= volCount[i];
			gradU[2][i] /= volCount[i];
		}
	}

	public double computeL2Norm(double[] u) {
		double result = 0.0;
		for (int i = 0; i < mesh.getNt(); i++) {
			int[] ele = mesh.getElement(i);
			mesh.getCoordPerElement(ele, coordRegion);
			double volume = CG.volume3DT4(coordRegion);
			double[] uu = MVO.getSlice(u, ele);
			result += volume
					* (uu[0] * uu[0] + uu[1] * uu[1] + uu[2] * uu[2] + uu[3] * uu[3] +
					   uu[0] * uu[1] + uu[0] * uu[2] + uu[0] * uu[3] +
					   uu[1] * uu[2] + uu[1] * uu[3] +
					   uu[2] * uu[3])
					/ 60;
		}
		return Math.sqrt(result);
	}

	public double computeH1SemiNorm(double[] u) {
		double result = 0.0;
		double[][] grad = new double[4][3];
		double dxU, dyU, dzU;
		for (int i = 0; i < mesh.getNt(); i++) {
			int[] ele = mesh.getElement(i);
			mesh.getCoordPerElement(ele, coordRegion);
			double volume = CG.volume3DT4(coordRegion);
			CG.gradient3DT4(coordRegion, grad);
			double[] uu = MVO.getSlice(u, ele);
			double[] dx = MVO.col(grad, 0);
			double[] dy = MVO.col(grad, 1);
			double[] dz = MVO.col(grad, 2);
			dxU = MVO.dot_product(dx, uu);
			dyU = MVO.dot_product(dy, uu);
			dzU = MVO.dot_product(dz, uu);
			result = result + volume * (dxU * dxU + dyU * dyU + dzU * dzU);
		}
		return Math.sqrt(result);
	}

	public double[] getNodeInfoFromElementInfo(double[] eleInfo) {
		int len = eleInfo.length;
		ParamCheck.checkEqual(len, this.mesh.getNt());
		double[][] coord = new double[4][3];
		double[][] grad = new double[4][3];
		double[] volCount = new double[len];
		double[] eleU = new double[4];
		double[] nodeInfo = new double[this.mesh.getNv()];
		for(int i = 0, n = mesh.getNt(); i < n; i++) {
			int[] ele = mesh.getElement(i);
			mesh.getCoordPerElement(ele, coord);
			double volume = CG.volume3DT4(coord);
			for(int j = 0; j < 4; j++) {
				nodeInfo[ele[j]] += eleInfo[i] * volume;
				volCount[ele[j]] += volume;
			}
		}

		for(int i = 0, n = mesh.getNv(); i < n; i++) {
			nodeInfo[i] /= volCount[i];
		}
		return nodeInfo;
	}
}
